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ABSTRACT 

We study and elucidate the mechanism of spiral density wave excitation in a differ- 
entially rotating flow with turbulence which could result from the magneto-rotational 
instability. We formulate a set of wave equations with sources that are only non-zero in 
the presence of turbulent fluctuations. We solve these in a shearing box domain, sub- 
ject to the boundary conditions of periodicity in shearing coordinates, using a WKBJ 
method. It is found that, for a particular azimuthal wave length, the wave excitation 
occurs through a sequence of regularly spaced swings during which the wave changes 
from leading to trailing form. This is a generic process that is expected to occur in 
shearing discs with turbulence. Trailing waves of equal amplitude propagating in op- 
posite directions are produced, both of which produce an outward angular momentum 
flux that we give expressions for as functions of the disc parameters and azimuthal 
wave length. 

By solving the wave amplitude equations numerically we justify the WKBJ ap- 
proach for a Keplerian rotation law for all parameter regimes of interest. In order to 
quantify the wave excitation completely the important wave source terms need to be 
specified. Assuming conditions of weak nonlinearity, these can be identified and are 
associated with a quantity related to the potential vorticity, being the only survivors 
in the linear regime. Under the additional assumption that the source has a flat power 
spectrum at long azimuthal wave lengths, the optimal azimuthal wave length pro- 
duced is found to be determined solely by the WKBJ response and is estimated to be 
2TrH, with H being the nominal disc scale height. In a following paper by Heinemann 
& Papaloizou, we perform direct three dimensional simulations and compare results 
manifesting the wave excitation process and its source with the assumptions made and 
the theory developed here in detail, finding excellent agreement. 
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1 INTRODUCTION 

Accretion discs are ubiquitous in astrophysics, occurring in 
close binary systems, active galactic nuclei and around pro- 



tostars (see e.g. Papaloizou & Lin 1995 Lin & Papaloizou 
|1996[ for reviews). Ever since their importance was first re- 
alized it has been clear that some form of turbulence is nec- 
essary to provide the anomalous angular momentum trans- 
port implied by observed luminosities and inferred accretion 
rates. This has usually been parametrised using the |Shakura| 
& Syunyaev ( 1973 1 a-parametrisation. 



The most likely source of turbulence is through the 
magneto-rotational instability (MRI) (see [Balbus &: Haw-] 
|ley|199l"| |1998[ ). Both local and global simulations that dis- 
play sustained MRI turbulence have been performed with 
and without net flux. In all cases proliflc spiral density (SD) 
wave excitation has been noted (e.g. jGardiner fc Stone|2005| 



in the local case and Armitage 1998 in the global case). 
These waves may be crucial for explaining various phenom- 
ena in accretion disk systems. In the context of protoplane- 
tary disks for instance, it has recently been suggested that 
stochastic gravitational forces derived from density varia- 
tions due to SD waves may play an important role in driv- 



ing the migration of low mass protoplanets (Nelson & Pa- 



|paloizou|2004| |Nelson|2005[ ) . This possibly remains a viable 
mechanism even in so-called dead zones where the ioniza- 
tion fraction is too low for MHD turbulence to occur, see 
the recent simulations by |Oishi, Mac Low fc Menou||2007| 
In general, SD waves may lead to significantly enhanced an- 
gular momentum transport in magnetically inactive regions 
of accretion disks. It is therefore important to gain an under- 
standing of the processes leading to the excitation of spiral 
density waves, how generic the phenomenon is, and how the 
wave amplitudes scale with physical parameters. 
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In order to do this, we assume weakly nonlinear condi- 
tions, under which the important source terms for exciting 
the SD waves are expected to be proportional to what we 
call the pseudo potential vorticity (PPV), which is equal to 
the potential vorticity to linear order but differs from it in 
the nonlinear regime. The potentially important role that 
PPV plays for the excitation of SD waves in rotating shear 
flows has been recognised in earlier work (see e.g. ,Chagel-, 
ishvili et al.||1997 l. These authors found by numerically in- 



tegrating the linearised equations of motion of compressible, 
plane Couette flow that there exists a linear mode coupling 
between specifled vortical perturbations and (free) SD waves 
that leads to efficient excitation of the latter (see also |Bodo| 
|et al.|[2b05| for a similar study in the context of accretion 
disks). However, the possibility of the generation of vortical 
perturbations through the action of turbulence has not yet 
been fully assessed. 

It is the purpose of this paper to develop a mathemat- 
ically rigorous theory of the excitation of SD waves within 
the WKBJ framework. Our approach is inspired by the work 



of Vanneste & Yavneh ( 2004 1 who developed an analogous 
theory for small amplitude inertia-gravity waves in a local, 
quasi-geostrophic model of the Earth's atmosphere. Apart 
from illuminating the excitation process of SD waves in ro- 
tating shear ffows by putting it on a firm mathematical ba- 
sis, the theory also enables us to calculate the amplitude 
of the excited waves and the associated angular momentum 
transport explicitly. 

In a following paper ( [Heinemann fc Papaloizou||2009[ 
paper II) we perform numerical simulations which to study 
SD wave excitation in MRI driven turbulence (see also |Gar^ 
[diner fc Stone||2005[ [Shen et al.]|2006| ) and make detailed 
comparisons to the WKBJ theory presented here. In this 
scenario, magnetically dominated turbulent stresses cause 
vortical perturbations which then in turn lead to the exci- 
tation of the observed SD waves. The role that turbulent 
stresses play is thus indirect. In this aspect the excitation 
process differs from that considered in Lighthill's theory of 
aerodynamic noise generation | |Lighthill[[l952| ). 

The plan of this paper is as follows: In section [2| we 
describe the shearing box model, giving the basic equations 
and defining the background shear flow in which the hydro- 
magnetic turbulence responsible for the SD wave excitation, 
is generated. In section[3|we derive equations describing the 
excitation of SD waves. These take the form of linear wave 
equations with both linear and nonlinear source terms that 
are determined by the turbulence. We focus on waves that 
are nearly independent of the vertical coordinate which have 
been found to dominate in simulations carried out in paper 
II and which can be dealt with using a vertical averaging 
procedure. We formulate the law of conservation of angular 
momentum for linear, non-dissipative SD waves and derive 
a convenient expression for the average radial flux which we 
later use to estimate the angular momentum flux arising due 
to the excited waves. 

In section |4| we go on to develop the WKBJ theory of 
wave excitation. A Fourier analysis is carried out enabling 
each azimuthal wave number ky to be considered separately. 
The WKBJ theory applies to a shearing box for which the 
boundary conditions are the imposition of periodicity in 
shearing coordinates and involves a sequence of excitations 
uniformly spaced and localized in time, during which the 



wave swings from leading to trailing. The wave amplitude 
and wave action produced in a swing are calculated from a 
WKBJ formalism involving the evaluation of integrals along 
anti-Stokes lines. 

We compare the results derived from asymptotic theory 
to results obtained by numerically integrating the ordinary 
differential equations describing the evolution of the appro- 
priate Fourier amplitude. We find excellent agreement be- 
tween these approaches. This agreement persists under all 
conditions of interest. This is in spite of the fact that asymp- 
totic theory formally requires a parameter depending on the 
azimuthal wave number to be small. Finally we discuss our 
results in Section |5l 



2 THE SHEARING BOX MODEL 
2.1 Basic set up and equations 

We consider a conducting gas in the shearing box approxi- 



mation ( Goldreich fc Lynden-BeTl[[1965 1. A Cartesian coor- 
dinate system {x, y, z) with origin at the centre of the box is 
adopted. The system rotates with angular velocity $7 = Q,ez, 
with being the unit vector in the z-direction. This coin- 
cides with the angular velocity of the centre of the box, taken 
to be in a circular orbit. In the Keplerian case this is about 
a central point mass. The lengths of the sides of the box in 
the three coordinate directions are {Lx, Ly,Lz) and vertical 
stratification is neglected. 

The basic equations are those of MHD for an isothermal 
gas, i.e. the continuity equation 



(p«) = 0, 



the momentum equation 



d{pv) 
dt 



= -c^Vp -inx pv- pV$ + V ■ T' + V • {2puS) 



and the induction equation 

dB 

dt 



V X (u X B - ?7V X B) 



where p is the density, the velocity is w = {vx,Vy,Vz), the 
isothermal sound speed is c, the magnetic field is i? = 
{Bx,By, Bz), the nonlinear stress tensor T' has components 



T' 



BiBj — SijB /2 — pviVj, 



and S is the traceless rate-of-strain tensor whose compo- 
nents are given by 



S^ 



2 \dxi 



dxi 



j5,wV ■ V 



The kinematic viscosity is v, the resistivity is rj, and the 
combined gravitational and centrifugal potential is given by 



o2 2 
-qil X , 



where for a Keplerian flow the constant q = 3/2. 

The isothermal MHD equations admit the definition of 
a characteristic length scale H = c/Q, which we will refer 
to as the nominal disc scale height even though we have 
neglected vertical stratification. 
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2.2 Equations for deviations from the steady state 
and boundary conditions 

The background state is taken to have uniform density po, 
zero magnetic field and a linear shear corresponding to the 
velocity 



-qClxBy 



As we are interested in wave propagation we work in 
terms of velocity deviations from the background shear, 
u = V — vo, which we use to define the linear momentum 
density per unit volume p = pu. In terms of p and p the 
governing equations now read 

Vp + V-p^O, (1) 

Vp = -c^p - 2n X p + q^p^Ey + V • T + V • {2puS), (2) 

VB = V X [ux B -rfj X B)- qflB^By (3) 
where the differential operator 

(4) 



accounts for advection by the linear shear. Note that here, 
the nonlinear stress tensor 



Tij = BiBj — SijB /2 — pUiUj 



(5) 



only contains products of deviations from the background 
state. 

We consider equations ([TJ to ([3]l to be subject to pe- 
riodic boundary conditions in 'Lagrangian' (or 'shearing') 
coordinates given by 



x' = X, y' ~ y + q^ltx, z = z, and = t, 



(6) 



transformation to which removes the explicit a;-dependence 
contained in ([4]| - albeit at the expense of an explicit 
time dependenc^ Re-expressed in 'Eulerian' coordinates 
(x, y, z, t) the radial boundary condition for any fluid vari- 
able / then reads 

f{x + L^,y- qfltL^, z, t) = /(x, y, z, t) 

while the azimuthal and vertical boundary conditions are 
simply 



and 



respectively. 



f{x, y + Ly,z,t) ^ f{x, y, z, t) 
f{x,y,z + L^,t) ^ f{x,y,z,t), 



3 WAVE EQUATIONS WITH SOURCES 

In order to proceed we develop equations for the deviation 
of the state variables from their background state in order to 
obtain wave equations with sources. These can then be used 
to study the excitation of SD waves explicitly. It is known 
that SD waves can propagate in a strictly isothermal box 



with no dependence on the vertical coordinate (see Fromang 



^ Witliout loss of generality we have assumed tfiat the two coor- 
dinate systems coincide at t = 0. 



|fc Papaloizou|2007[ ) and we have found that such waves are 
the ones predominantly excited in our simulations presented 
in paper II. We therefore vertically average equations ([T]) and 
([2| , which then become equations for the vertical averages of 
the state variables such as p and p, in order to describe such 
waves. When periodic boundary conditions in z are adopted 
this can be done without approximation and it is equivalent 
to adopting kz — when Fourier transforms are considered. 
Proceeding in this way, we denote the vertical average of a 
quantity / by use of angle brackets as (f) ^■ 

To further simplify the analysis we consider the inviscid 
limit of the shearing box equations, which then become 

V{p)^+d^{p,.)^ + dy{py)^=0 (7a) 

V (p,)^ + c^d^ (p)^ - 2n (py)^ = (01,)^ (7b) 



1^ (Py). + cdy (p), + (2 - q)Q (p,)^ = (my)^ 



(7c) 



where we introduced the short-hand — V-T. At this point 
we note that in the zero net flux case considered here, the 
magnetic fleld enters the momentum equation only through 
the nonlinear stress tensor ([5| so that it will not affect the 
description of linear SD waves. 

Acting on ( |7c| | with V and rearranging terms yields 

{V'-C'\7' + K'){py)^ = 

-c^d.{0z + {q^2)n{m.)^+V{yiy)^ (8) 

where = 2(2— g)r2^ is the square of the epicyclic frequency 
and 

C = dxPy - dyPx + {q- '2)Qp, 

which we shall call the pseudo potential vorticity (PPV). To 
linear order, the variation of PPV is equal to the variation 
of potential vorticity (PV), 



xUy — dyUx + (2 — q)0. 



(see [Johnson fc Gammie||2005[ ). Written out explicitly, we 
have 



— (^0 ~ PoiQ ^ Qo) + norflinear terms. 



where 



Co = ((? - 2)flpo and Qo = 



(2 - q)Q, 



pa 



are the steady state background values of PPV and PV, 
respectively. For disturbances with fe^ = and given a 
barotropic equation of state, PV is an exactly conserved 
quantity whereas PPV varies due to nonlinear stresses, 

T^{0.=dx{^y)^~dy{^x)^. (9) 

We can form wave equations similar to ([8| for p and p^ ■ 



Letting V act on ( 7a I and ( 7b I yields 



(©2 _ ^2^2 ^ ^2-| ^^^^ ^ 2gf^a^ (p,)^ = 



2n{c),-a.(ai.),-a,(0T,), (loa) 



and 



c^dy{0^ + 2n{^y)^ + v{<nx)^. (lOb) 
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3.1 Decomposition into shearing waves 

In the periodic shearing sheet we may expand all fluid vari- 
ables in a series of plane wave solutions 



exp 



(^ik'^x' + ik'yy' 



in the Lagrangian coordinate frame ([sjl . Here, the radial and 
azimuthal wave numbers 

, I 2-Kny 

— — — and ky = — — - with nx,ny £ £. 

In terms of Eulerian coordinates, the plane wave expansion 
for any (vertically averaged) fluid variable (/)^ reads 

(a;, y, i) = ^ f{t) exp [ifc^(i)i + ifeyj/] 

with a time dependent radial wave number 

kx{t) = k'x + qntk'y. (11) 

and constant azimuthal wave number ky = ky. 

When viewed from the Eulerian coordinate frame, the 
radial wave number of non-axisymmetric disturbances (for 
which ky 7^ 0) changes linearly in time due to advection 
by the linear shear, leading to the notion of sheared distur- 
bances as originally discussed by Kelvin ( |Thomson|| 1887[ ) . 
In the shearing sheet, non-axisymmetric plane waves are 
therefore often referred to as shearing waves. In an astro- 
physical context, the usefulness of the concept of sheared 
disturbances for understanding accretion disk phenomena 



was first realised by Goldreich & Lynden-Bell (19651. 

It is customary to classify shearing waves according to 
whether they are leading, i.e. kx{t)/ky < 0, or trailing, i.e. 
kx{t)/ky > 0. Because the time dependence of the radial 
wave number (111 is such that kx{t)/ky always increases 



monotonically provided that gil > 0, every leading wave will 
eventually become trailing as time progresses. The change 
from leading to trailing is referred to as 'swing' and occurs 
when kx(t) — 0. Different shearing waves swing from leading 
to trailing at different times, successive swings being sepa- 
rated, for a given fey, by a flxed time interval 



St, = 



(12) 



where Torb ~ is the orbital period. 

We note that because we are dealing with the Fourier 
transforms of real quantities, each Fourier coefficient be- 
comes equal to its complex conjugate under reflection of 
the wave number such that fc — > — k. This means that we 
may without loss of generality consider only ky > 0, and 
then multiply by a factor of two after taking the real part 
of the transforms to obtain physical quantities, which then 
accounts for ky < 0. Thus from now on we consider only 

ky > 0. 

It is now straightforward to write down the SD wave 
equations ([8| and ( |10[ | in Fourier space. Substituting 



2? d/dt, dx ikx{t), and dy 



iky 



we readily obtain 



d^ 
dt2 



p + 2qQ,\kyPx — 
- 2Q.C, - ikx{t)Vix - iky^y, 



(13a) 



+ [k'^{t)c^ + K^] px + 2qQc^ikyp = 



dt^ 



and 

d^Py 



+ [k\t)c' + K^] py = -c^ikxC + {q- 2)nmx + 

(13c) 



c\ky( + 2Qmy + ^, (13b) 



dmy 



3.2 Angular momentum flux 

Conservation of angular momentum for linear waves in the 
shearing box follows from invariance of the system under 
translations along the azimuthal or y-direction. Here we note 
that this actually yields a momentum flux that can be con- 
verted into an angular momentum flux by multiplying by 
the radius of the centre of the box. As the latter quantity 
does not play any role in the box dynamics we can conve- 
niently set it to be unity making the momentum and angular 
momentum fluxes equivalent. We introduce the Lagrangian 
displacement which we define through 



2?^ = p/po- q^^xBy 



(14) 



In terms of ^, the density deviations from the background 
state are given by 

Sp/po + V ■ I = 
and the linearised equations of motion become 

V^S. = cVV ■ I - 2n X 2?^ + 2qn'^£,xex- 
These also follow from the requirement that the action 



S = 



Jc{i,Vt\7-i) d'xdt, 



(15) 



S ^ J Cd'xdt, 
with the Lagrangian density given by 

c=^[vi-vi-c\\/- if + 2{n xi)-vi + 2qn'e.] , 

and the integral being taken over the box and between two 
arbitrary points in time, be stationary with respect to arbi- 
trary variations of the Lagrangian displacement. The angu- 
lar momentum conservation law follows from the invariance 
of the action ( |15[ ) under infinitesimal translations in the y- 
direction. The resulting form of Noether's theorem yields 



V 



dC 



dC 



= 0. 



We thus define the angular momentum density 
dC 



A 



dm) 

and the angular momentum flux 
dC 



dyi = ^POm + flXi)-dyi 



(16) 



9(V-0 ' 
which enables us to write 

VA + V ■ F ^ 



dyi + Cey = Spdyi + Cey, (17) 



(18) 



Note the minus sign in equations ( |16[ ) and ( |17[ ) which can be 
determined from considering the action of an external force 
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(see also |Ryu fc Goodmcin||1992[ ) . The angular momentum 
conservation law ( 18 1 may be averaged over y and z to yield 



dt{A)y,+d:,{F^)y, =0. 



The angular momentum flux as defined in ( 18 1 is, at 
times, inconvenient to work with because it involves the La- 
grangian displacement ^. We can derive a related wave ac- 
tion where the radial fiux only depends on mass density and 
momentum density as follows. From \7c\ and ( |14[ | we have 
to linear order 

Py + C^dyW 

(q - 2)npo ' 

Inserting ( |19[ ) into ( |17| l yields after some straightforward 
algebra 



where VW = 5p. 



(19) 



(2 - q)npo 

SpdyPy+dy(^C^5pdyW^ 



V 



dyW)^ I 



The second and the third term in curly brackets can be 
absorbed in the j/-component of the angular momentum fiux 
and in the angular momentum density, respectively, giving 
rise to a new wave action conservation law, 

VA' + V • F' = 0, 

where the modified radial angular momentum flux 

, ^ c^SpdyPy 
" (2 - q)npo ■ 

has the desired property that it does not involve the La- 
grangian displacement. Furthermore, this new flux is equal 
to the radial component of the angular momentum flux ( |17[ ) 
after averaging over y, z, and t. But it should be noted that 
it in order to establish this equality it has been assumed 
that no external forces act in the domain. 

We will define a suitable temporal averaging procedure 
further below in Section |4.11[ At this point we note that 
when evaluated for a single pair of (complex conjugate) 
shearing waves, and averaged over both y and z, the two 
equivalent expressions for the radial angular momentum fiux 
become 



and 



(2 - q)flpo 



Im(p*p 



(20) 



(21) 



Here, without loss of generality, we have adopted fcj, > as 
described in Section |3JJ 



4 WKBJ THEORY OF WAVE EXCITATION 

In this section we derive a WKBJ theory of the wave exci- 
tation that occurs during a swing cycle and derive an ex- 
pression for the wave action produced. We go on to com- 
pare this theory in detail with the results of numerical inte- 
grations of the ordinary differential equations governing the 
time dependent evolution of the Fourier transforms of the 
wave amplitudes. Excellent agreement is obtained. In paper 
II we compare results obtained from the WKBJ theory with 
those obtained from MRI simulations. 



4.1 The nature of the source terms 

We first need to establish which of the source terms on the 



right hand sides of ( 13 1 are primarily responsible for wave ex- 



citation. Inspection of these equations shows that the source 
terms are of two kinds. The first kind is proportional to the 
transform of PPV and the second kind is proportional to 
the nonlinear stress tensor. Only terms of the first kind re- 
main in the linear regime. Thus under conditions of weak 
nonlinearity, we would expect them to dominate. An analy- 
sis of the relative contributions found in direct simulations 
given in paper II shows that the contribution of the pseudo 
potential vorticity related terms is the more important by 
an order of magnitude confirming the above idea. It is also 
shown that the strength of the wave excitation phenomenon 
for a swinging wave is directly correlated with the ampli- 
tude of the pseudo potential vorticity transform at the time 
of the transition from leading to trailing. Here we reiterate 
that although the generation of PPV itself is driven by the 
nonlinear stresses, see (|9|, the linear source terms involving 
PPV survive at linear order if there is a build up over time 
under conditions of weak nonlinearity. 

Following on from the above discussion, from now on we 
retain only source terms that depend on the pseudo potential 
vorticity. These are proportional to which we recall is a 
Fourier amplitude which is given by 



1 



LxLy JQ JQ 



From the above expression we note that if {Q ^ {x, y, t) is an 
ultimately smooth function, the Riemann Lebesgue lemma 
allows us to infer that the source is negligible as t ^ ±oo and 
so is expected to peak when the wave swings from leading to 
trailing. Thus the wave amplitude excitation process should 
also be localized around this time. 

At this point we recall that in a shearing box the wave 
excitation process for a fixed ky appears as a succession of 



swings, which as seen from (111 are separated by the time 
interval ( |12[ ). For the longest possible wave length in the y 
direction given by the box size this is Ly / {2-nqLx) orbital 
periods. This is independent of the box size as long as the 
aspect ratio is fixed. Although this time interval formally 
decreases as increases, because of the periodic symmetry 
associated with the shearing box, we expect results to be 
independent of once this is larger than the radial corre- 
lation length associated with the turbulence, expected to be 
~ H. Accordingly we might expect the phenomenon to take 
a similar form in global simulations (e.g. |Nelson|2005| . 



4.2 Reduction to three uncoupled second order 
oscillator equations 

As motivated above we will now neglect any nonlinearities in 
the problem. After dropping the nonlinear source terms ap- 
pearing in the SD wave equations Q and ( 10 1 the evolution 
of p, px, and Py in Fourier space is governed by 



^ + {k'c +k'^p + 2qnikyPx = -2nC, (22a) 



^ + [k^c + K^'^px + 2q^c^ikyp = cikyC, (22b) 
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and 



4.3 Slowly varying solutions 



(23) 



We note that in the absence of nonlinearities the pseudo 
potential vorticity, 



is conserved exactly, i.e. 



dt 



= 0. 



(24) 



(25) 



The wave equation for py, given by ( |23|l, therefore decou- 
ples from those for p and Px, given by (221. The latter two 



equations may be decoupled from each other as well by in- 
troducing the pair 



P± ^ Px± pC, 

for which the linearised wave equations read 



(26) 



d p± 



+ (k 



-'■^c^ + k"^ ±2qQikyCjp± = c(ikyC^^ 2Q^C- (27) 

We reiterate that the above equations describe the excitation 
of density waves and accordingly it may be confirmed that 
they are absent in the incompressible limit for which c ^ oo 
while maintaining time derivatives finite. In this limit, it is 



readily verified from ( 22 1 and ( 23 1 that the excited or forced 



flow considered below simply becomes the incompressible 
flow associated with a vorticity distribution that is slowly 
varying through the action of nonlinear MHD forces. We 
note the pseudo potential vorticity becomes the vorticity in 
that limit. 



We can simplify ( 23 \ and ( 27 \ further by introducing 



the dimensionless time variable 

kxC 

and the dimensionless parameter 

q^lkyC 



in terms of which we have 



2 d^Py 

^ dT2 



ikxC 



(28) 



(29) 



(30a) 



and 



d P± 
dr2 



1 ± 2ie)p± = c 



ikyC =F 2n\ 
kyC^ + k'^ J 



C. (30b) 



We remark that the homogeneous form of these equa- 
tions with C = can be solved in terms of Parabolic Cylin- 
der functions (e.g. |Narayan et al.|1987| |. However, we did not 
find this feature to be useful in the context of this paper. 
Rather we use the fact that the form of the inhomogeneous 
equations suggests an asymptotic expansion in e. Formally, 
such an expansion is valid if e ^ 1. From ( |29| l we can see 
that this will be the case both in the high and the low az- 
imuthal wave number limit given by kyC 2> k and kyC <^ k, 
respectively. 

In the following, we will derive asymptotic solutions to 
| |30[ ) based on the smallness of e. Ifowever, we will demon- 
strate that these approximative solutions show excellent 
agreement with the exact solution obtained from direct nu- 
merical integration even in the worst possible case oikyC = hi 



We consider equations ( 30 1 for large jr| , i.e. in the high radial 
wave number limit. In this limit the second time derivatives 
are significant only for solutions that vary rapidly in r. Such 
solutions are expected when sources are absent and corre- 
spond to high frequency oscillations. This suggests that they 
can be dropped for solutions that vary slowly with r. As a 
first approximation, we drop the double time derivative in 
( 30 1 to obtain 



and 



Py 



p± = c 



ikxC 



ikyC =F 20, 



c 



fc^c^ + k2 ± 2qQ.\kyC 



(31a) 



(31b) 



where we have used ( 28 1 and ( 29 1 . These approximative solu- 



tions to inhomogeneous problem are the leading order terms 
of an asymptotic series expansion in ascending powers of e. 
The series diverges near the time of the swing from leading 
to trailing, i.e. near r = 0, when the double time derivative 
in ( |30| l becomes significant and oscillatory solutions to the 
homogeneous equation must be taken into account. 

The occurrence of such oscillatory solutions may eas- 
ily be demonstrated by direct numerical integration of ( |30[ ), 
which is just a set of linear ordinary differential equations. 
We start the integration in the far leading phase, i.e. at 
r = To with To large and negative. In this limit the balanced 



solutions (31 1 hold and may be used as initial conditions. 



We note that in doing so we have to be careful not to vio- 
late PPV conservation which we know to be exact in linear 
theory, see (25 1. This problem arises on account of the ad- 



ditional time derivative taken to obtain ([8| and 1 10 1 from 

Formally, the balanced solutions ( |31| l are reconcilable 
with PPV conservation only in the limit r — + ±oo. We thus 
introduce an error if we use these solutions as initial con- 
ditions at some finite tq and we have to make sure that 
To is sufficiently large so that this error is small. In order 
to be able to quantify this error during the course of the 
integration, we express the PPV C, in terms of p and px, 



for which e attains its maximum value, e" 



qQ./2K. 



see (241 together with (26 1, and solve (30 1 as a system of 
three coupled ordinary differential equations. (Alternatively, 
equations ([7| could be solved directly). PPV conservation 
is not guaranteed in this case, but we find empirically that 
the numerical integration conserves PPV arbitrarily well de- 
pending on how large ro (and thus the error introduced by 
using the balanced solutions as initial conditions) is. 

Bearing these general remarks in mind, we now discuss 
a specific numerical solution to ( |30[ ). For this purpose we 
consider Keplerian shear, i.e. q = 3/2. Because we would 
like to determine empirically how well asymptotic theory 
works when we are far away from the asymptotic limit e ^ 1 
we take the worst possible case, i.e. a shearing wave with 
ky = k/c= 1/H so that e — qQ./2Hi = 3/4. For the determi- 
nation of the initial conditions from the balanced solutions 
we assume, without loss of generality, that Cj= ^po in (31 1. 

We start the integration at r = — lOOPlThe evolution 



^ At this point the relative error as far as PPV conservation is 
concerned is ~ 10~*. We find that the relative error never exceeds 
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Re p/po 




Inipy/ pQC 




Numerical 
Balanced 



Figure 1. Comparison between the numerical solution (black) 
to the linearised wave equations (jSOjl with the corresponding bal- 
anced solutions (blue). Here we have used | |26| to compute p and 
from p±. The parameters in this example are q = 3/2, e = 3/4, 
and = f2po. 



of the Fourier amplitudes near the swing from leading to 
trailing as a function of the radial wave number is shown in 
Fig. [l] Before the swing from leading to trailing the numer- 
ical solution closely follows the balanced solution up until 
T « where a sudden transition to oscillatory behaviour 
occurs. In the trailing phase, the numerical solution oscil- 
lates with a relatively large but rather slowly evolving am- 
plitude around the balanced solution. Such oscillatory be- 
haviour cannot be captured by a perturbation series expan- 
sion. In order to describe it we have to resort to singular 
perturbation theory to be discussed in the next section. 

4.4 WKBJ solution of the generic oscillator 
equation 

To study the excitation of a general Fourier mode we seek 
a solution to the forced harmonic oscillator equation of the 
general type 



e + (r + a )j/(r) = /(r) 



(32) 



10 ^ during the entire course of the integration from r = —100 
to T = 10. 



with e ^ 1 and where a is a complex number with 
\a\ ~ Oil) and |pha| < 7r/4, 



(33) 



the latter being true for a — \/l 4- 2ie with ^ e < oo. 

Throughout the following analysis make w = z^^" 
single-valued by taking the branch cut along the negative 
real axis and always take the principal branch, thus 



phw 



ph z 



(34) 



We note that equations ( |30a[ ) and ([30b| are special 
cases of the above general form. We further also assume 
that /(r) is a slowly varying function such that A/ ~ 0(1) 
for At ~ 0{1) as e 0. Physically this means that the 
pseudo potential vorticity transform should vary at a slow 
rate compared to the wave oscillation frequency at around 
the time of the swing. This is expected under conditions of 
weak nonlinearity as discussed in section [44] and simulation 
results presented in paper II indicate that this is indeed the 



4.5 Outer and balanced solutions 



A solution of equation ( 32 1 in ascending powers of e is read- 



ily found from regular perturbation theory, and yields to 
lowest order 



y{r) = 



fir) 
-2 + 



+ 0{e 



(35) 



which is, given (331, well defined on the entire real axis. 



However, this series, giving rise to what is described as the 
balanced solution, is only asymptotic. At all orders in e^, 
the series will fail to represent rapidly changing oscillatory 
contributions to i/(r) for which the double derivative term 
in (321 is significant. These vanish more rapidly than any 
power of e as e — > 0. 

In order to obtain these contributions, that are in fact 
associated with the excited SD waves, we have to resort to 
singular perturbation theory. We thus seek a solution of the 
asymptotic form 

f(r) ^ , ^„ 



yir) 



yij) 



fir) 



+ 



A+e' 



i#(r)/E 



+ 



A_e+' 



(r2-Ha2)i/4 (^2,^^2)1/4' 



(36) 
(37) 



which is the sum of the leading order solution to the in- 
homogeneous problem obtained from regular perturbation 
theory, and a standard WKBJ solution to the homogeneous 
problem. Here the WKBJ phase is given by 



$(r) = / \/cr2 + a2 da 



T\/r2~+~o2 _|_ 



+ + a2 



(38) 



4.6 Matching on anti-Stokes lines 

In order to determine the WKBJ (or wave) amplitudes A+ 
and A-, we analytically continue the WKBJ solution ( |37[ ) 
into the complex r-plane and match it to inner solutions 
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valid in the immediate vicinity of the two complex WKBJ 
turning points at r, — ±ia. Thus in this section the ± alter- 
native refers to the turning point with positive and negative 
imaginary part respectively rather than the p± alternative 
of the previous section. 

We find it convenient to match inner and outer solutions 
along the so-called anti-Stokes lines defined by 

lm$(r) = Im<l>(r.). (39) 

Because the imaginary part of the WKBJ phase is constant 



ph a = 



ph a = tt/ 



on these lines, one of the WKBJ exponentials in (371 is 
maximally sub-dominant to the other one and may thus be 
ignored. 

From the first order Taylor expansion about a turning 
point we have 



2 I 2 

T + a 



2r.(r- r.) 



(40) 



and therefore from equation (|38|l we obtain 

n3/2 

r,) 



$(r) ^ $(r* 



2^2 1/2 



The condition ( 39 1 thus defines three anti-Stokes lines ema- 



nating from each turning point at angles given by 

5tv tt , 7r\ pha 



ph(r 



6 



IT TV 

^6' ±2 



For consistency we match on the anti-Stokes lines originating 
from each turning point that asymptotically approach the 
real axis. 

For Rer < Re(±ia) these are 



ph(r - 

and for Rer > Re(±ia) 
ph(r- 



5lT 



pho 

~3~ 



(41) 



(42) 



TV ph a 

Note that with the definition of the square root obtained 
from ( |34[ ), the WKBJ solution ( |37[ ) has branch cuts that 
leave the turning points at an angle 



ph(r 



ph a 



so we can always match along the anti-Stokes lines specified 
above without crossing branch cuts. This is illustrated in 
Fig.[2l 

On an anti-Stokes line one of the two WKBJ exponen- 
tials in ( |37[ ) will be maximally sub-dominant to the other 
in the limit e ^ 0, depending on the sign of the imaginary 
part of the WKBJ phase at the turning point from which it 
emanates. At these we have 

<I>(±ia) = ±ia^7r/4 

and thus when the conditions ( |33[ ) are satisfied 

Im$(±ia) § 0. 

This means in general that the maximally sub-dominant 
exponential has an amplitude smaller by a factor oc 
exp(— co/e), where Co is a constant of order unity when com- 
pared to the dominant one, which is very small for small 
e. Therefore it should be neglected. Dropping the maxi- 
mally sub-dominant WKBJ exponential, in the vicinity of 



1^ 

6 





1 

Jl 












r 

1 







Rer 




Rer 



Figure 2. Anti-Stokes lines (solid black) emanating from the 
two turning points for pha = (left panel) and pha = tv/S (right 
panel). The curly red curves indicate branch cuts. 



pho- = +2tv/3 



pha = —2tv/3 



1 

B 
-1 



























Re z 



Figure 3. Integration contours for the Airy-integral ( |47| l in the 
z plane. In each case the contour is deformed from the positive 
real axis to the two paths shown. These contours (black), which 
extend to infinity, correspond to the situation where we are on 
the anti-Stokes lines that radiate away from the turning points 
towards the positive real axis on which ph c = 2tv/3 at the upper 
turning point (left panel) and phtr = —2iv/3 at the lower turning 
point (right panel). 



the turning points our WKBJ solution is to leading order in 
r =F ia 



A±e 



TTa 1 4e 



2r.(r - r,) 



[2r.(r-r,)]i/4 

^12^2 
3 



exp 



Tie 



3/2 



We see that matching on the anti-Stokes line emanating from 
the upper (lower) turning point will only enable us to deter- 
mine the WKBJ amplitude A+ (^-). We will thus need to 
match on the anti-Stokes lines from both turning points in 
order to determine the full WKBJ solution. 



4.7 Inner solution 

To obtain the inner solutions we consider the governing 



equation (32 1 in the vicinity of the turning points. Using 
the first order Taylor expansion around a turning point ( |40[ ) 
this becomes 



e^'^+2r.(r- 



dr2 



Oy(r) = /(r.). 



(43) 
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We now defined rescaled variables througli 



g±2Ti/3 



(2r.) 



1/3^-2/3 



(44) 



j/(r) = e-^-'^^(2r*)-^/^e-^/^f(r,)0(<T) (45) 

The first of these indicates that the solutions we seek vary 
on a scale |r| ~ e^^^ <^ 1. This feature enables us to perform 
an asymptotic expansion valid for large \a\ and still remain 
in the vicinity of the turning points. In terms of the rescaled 



variables, equation ( 43 1 yields an inhomogeneous Airy-type 
equation 



dcr2 



(46) 



We remark that from equation ( 44 1 we deduce that 

1 n r T./ \ Stt pha 
ph(T = tor ph(r — r, ) = =1=^ 



which applies on the anti-Stokes lines given by (411 and 



Phcr = ±— for ph(r-r.)==F- — 

which applies on the anti-Stokes lines given by ( |42[ ). 



4.8 Integral representation for the inner solutions 



The solution to ( 46 1 can be written as an integral 
'dz with h{a, z) — — o 



z73. (47) 



We are interested in the large-cr asymptotic behaviour of 
this solution, which is applicable on each of the anti-Stokes 
lines emanating from both turning points. 

For the anti-Stokes lines corresponding to phcr = we 
integrate along the positive real axis and get an end-point 
contribution via Watson's lemma 

oo, ph (7 = 0) 



0(a) 



or with (Uil and Jia 



2t,{t - T,) 



which matches our outer solution ( |37[ ) to the left of the turn- 
ing points if there we set Aj^ = A- = consistent with 
the causality requirement of the lack of existence of excited 
waves for r — oo. 

For the anti-Stokes lines corresponding to 
ph(T = ±27r/3 we deform the integration contour from 
the positive real axis to a contour consisting of two separate 
line segments, writing 



/■oo 

Jo 



dz = 



dz + 



[ e'''"'"'dz. (48) 



The first integral on the right hand side is evaluated along 
a straight line that goes from the origin to oo e^'^'^'''''. Along 
this line, h{z) is purely real and negative, and we again get 
an end-point contribution via Watson's lemma. 



L 



pT2i,r/3 



dz ' 





(a. 



oo, phcr = ±27r/3). 
The curve r± along which the second integral on the right 



hand side of (481 is taken goes from ooe^^'^'''^ along a path 
of steepest descent through the saddle point Zs = ^^ia^^^ 
and from there to oo, see Fig. [3] In the limit cr ^ oo, most 
of the contribution to this integral will come from near the 
saddle point Za, at which h(zs) = ±(2/3)i(T^^^ and so 



dz 



^l/2g±(2/3)i<T3/2 



We thus have 



e=F-'/4crl/4 

> OO, pha = ±27r/3) . 

^l/2g±(2/3)icT3/2 

a ^ eTi'^/4ai/4 

> oo, pha = ±27r/3) . 



Using ( 44 \ and ( 45 I we may reexpress the above solution in 



terms of y and r — r* to obtain 



^±37ri/4 



2r,(T- 



+ 



(2n)3/4(T- r,)i/4 
1 2%/2_i/2 



exp 



Tie 



n3/2 



We match this solution to our outer solution ( |37[ ) in the 
neighbourhood of each turning point by equating the factors 
in front of the maximally dominant WKBJ exponentials. 
This then yields the amplitudes 



A± = ±i. 



/(±i«) ( 



TT \l/2 

I 



.2aeJ 

Using these, the full asymptotic solution takes the form 

1/2 „-ira^/4£ 



fir) 



1 
2i 



-I- 



(r2 -Ha2)i/4 
^_ e+i*(-)/^ _ /(+i„) e-'*(-)/'] . (49) 



4.9 Determination of the wavelike forms of Py, p±, 

Px, and p 

We are now in a position to obtain explicit solutions for py 
and p±. To bring the wave equation | |30a[ | for py into the 
form of the general oscillator equation (321 we set 



1 and fir) 



To obtain the wavelike part of the solution, i.e. the compo- 
nent proportional to the WKBJ exponential, which we will 
denote by a tilde, /(r) has to be evaluated at the complex 
turning point r,. Because PPV is conserved in linear theory 
this poses no difficulties and we immediately obtain 

-7r/4e 



Py 



iCsC 



^fe2c2 + k2 V e (r2-Hl)i/4 



i[<I.(r)/6], (50) 



where (s = C(0) denotes the PPV amplitude and the time 
of the swing and the WKBJ phase 

*(^) = I [rVr'^ + 1 + In + v/r2 + l)] . (51) 

Using the definitions of r and e, respectively given by (28 1 



and ( 29 1 , we can rewrite the solution for Py as 



Py — iHC,sA\/0./uj cos 



Jo 



LO dt' 



(52) 
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where the dimensionless ampUtude 



Rep/po 



A = 



ir/4e 



(53) 



and we have defined the SD wave frequency 

This procedure is readily repeated to obtain exphcit 
solutions for p±. In this case, comparison of the governing 
equation (|30b|) with (|32|| implies that we need to set 



Vl ± 2ie and /(r) = c 



c 



The solutions for p± determined from (|49| are thus 
'2Q^^ikyc\ 1 



p± = iCsC 



(l±2ie)i/4 

27r 



e (r2 + 1 ± 2ie) 
where the WKBJ phase 



T7isin(^$±/e), (54) 



IK 



t2 + 1 ± 2ie 



(1 ± 2ie)ln 



1 ± 2ie 



VI ± 2ie 

It is important to note here that $± are not real with the 
consequence that quantities such as e'*+''' have a power law 
as well as exponential dependence on r for large r. 

Again, we rewrite ( |54[ ) in more familiar terms using ( |28[ ) 
and ((291 to obtain 



(55) 



p± — i(sHA±\i^il/uj± sin tj± dt' 
where this time the amplitude is given by 



A± = 



(fc2c2+K2±2gnifcyc)i/4 



20, =F ifcyC 



ir/4e 



(56) 



and the wave frequency is 

Lj± = \/ LO^ i: 2qOikyC. 

Note that both the amplitude and the wave frequency are 
complex valued. 

Recalling that p-j- = p^ zt pc, it is now a simple matter 
to determine the wavelike contributions to p^ and p in the 
forms 



p — —(sO ^ Im 



LU+ dt' 



and 



Px = iCsHRe 



A+ \/Q,/uj+ sin uj-f- dt' 



(57a 



(57b) 



where we have used the fact that A+ = .41 and tj+ — ujI . 

We have seen in Section|4j3]that the balanced solutions 
| |31[ ) show good agreement with the exact numerical solution 
in the leading phase but fail to capture the oscillatory be- 
haviour in the trailing phase. Having derived the wave like 
WKBJ solutions of the homogeneous wave equations, i.e. 




Impy/poc 



1 1 

Numerical 
— Balanced + WKB 




Figure 4. Comparison between the numerical solution (black) 
to the linearised wave equations (jSOj with the full asymptotic, 
i.e. balanced plus WKBJ, solution (red). The parameters are the 
same as in Fig. [l] i.e. q = 3/2, e = 3/4, and ( = Qpo- In the 
lowermost panel, py was obtained from |58| , while in the panel 
immediately above it was obtained from 1521. 



( |52[ ) , ( |57a[ ) and ( |57b[ ) , we are now in a position to determine 
whether they correctly describe this oscillatory behaviour. 

In Fig. |4] we compare the exact numerical solution dis- 
cussed in Section [4. 3| with the full asymptotic solution, i.e. 
the sum of the balanced and the WKBJ solution, in the trail- 
ing phase (r > 0). There is excellent agreement between 
the full asymptotic solutions and the numerical solutions. 
Beyond t — 2 the asymptotic solutions are virtually indis- 
tinguishable from the numerical solutions. Remarkably, this 
is so even though for this example we have chosen a shear- 
ing wave with an intermediate azimuthal wave number of 
kyC = K for which the 'small' parameter e attains its maxi- 
mum value, e = 3/4 in the case of Keplerian shear considered 
here, and we are therefore as far away as possible from the 
asymptotic limit e <§; 1. 
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We note that in the case of py there is a small but not- 
icable discrepancy between the numerical solution and the 



full asymptotic solution obtained directly from (52 1. How- 
ever, we can derive an alternative expression for py from 
PPV conservation. Because the WKBJ solutions are solu- 
tions to the free wave equations they should carry no PPV, 
from which it follows that 

ikyPx + (2 — q)^p 



Py 



(58) 



This expression agrees with ( 52 1 in the limit e ^ 1 and we 
see from Fig. |4] that ( |58[ ) is in fact more accurate for e ~ 1. 

We comment that after we have reintroduced the spa- 
tial dependence by multiplying with exp(ife • x) and then 
taking the real part, these solutions are found to consist of 
two waves of equal amplitude travelling in opposite direc- 
tions. This is a natural outcome given the symmetries of 
the shearing box. However, for the same reason, both waves 
transport angular momentum in the same direction, i.e. out- 
ward if they are trailing, see Section |4.11| 



4.10 



Asymptotic behaviour of the WKBJ 
solutions 



The WKBJ solutions for p and p^ given by (57a I and (57b I, 



respectively, involve complex phases which disguises their 
large time asymptotic behaviour. To make this more appar- 
ent we note that for large times or equivalently r we have 



$4 



2 + 111(2^) 



(1 + 2ie) ln(l + 2ie) 



where $ is given by equation ( |51| l and is purely real. We see 
that the imaginary part of will result in an extra power 
of r when taking the exponential. This has the consequence 
that when the sine of the WKBJ phase in ( |57a[ ) and ( |57b[ ) is 
re-expressed in terms of exponentials only those with abso- 
lute values that increase with r need to be retained. In this 
case these are oc exp(—i$-|- /e) and their asymptotic form is 
given by 



exp^— i$+/e^ ~ 2r exp^— i<3>/e 
1 

VI + 2ie 
from which it follows that 



exp 



ln(l + 2ie) 



4ie 



{!l''^''')^'0 '''''' {-'ly') 



^4 

A, 



where we have defined 



: exp 



ln(l + 2ie) 
4i(^ 



(59) 



yJk'yC^ -\- 1<? -\- 2qQikyC 

Using the above relations we can readily find the large r 
asymptotic forms of p and Pa, given by 



-Csfi^^lS+lv^fi cos udt' -phB^ 



(60a) 



and 



p^ ~ iCs-ff|S+|\/^I77^^ sin udt' - phS+^ . (60b) 
We thus see that because the WKBJ amplitudes for p and 



Px are complex valued, see (561, the envelope of the oscilla- 



tion grows in time and there is non-trivial phase shift with 
respect to py. 



4.11 The angular momentum flux 



In section [3^ we obtained two equivalent expressions for the 
radial angular momentum flux which uses the radial La- 
grangian displacement and which uses the y-component 
of the momentum density. We determined the angular mo- 
mentum flux for a single pair of (complex conjugate) shear- 
ing waves. Here, we are interested in the angular momentum 
flux associated with the excited waves, i.e. with the WKBJ 



solutions. In (20 1 and (211 we therefore replace p, Py, by 
P, Py, ^x, respectively, and obtain 



and 



{Fx)yz 



(2 - q)npo 



Im 



PyP 



(61) 



(62) 



where fc^ > is understood. 

We can use the WKBJ solutions just obtained to deter- 
mine these angular momentum fluxes. In order to calculate 
( 61 1 we need an explicit expression for the radial Lagrangian 
displacement which in Fourier space is given by 

ix = — Px dt. 
Po J 

the leading order WKBJ solution in the limit of large r can 
be calculated from (|60b| directly and is given by 



cos I / Ludt' — ph 



Using this result together with ( 60a I we readily obtain 

kyc\U^H^\B+\^ 



lim (Fx),, 

With the help of ([56| and 
form 

lim JF^ 



Q,po 

this may be expressed in the 



-7r/2E 



qpo 



m 



{2qnkycY 



3/4 ' 

(63) 

where the overline denotes an average over one oscillation 
period, and we have deflned 

tan~^(2e) 



1 



2e 



Alternatively, we may use equation ( 62 1 to determine 



the wave action. This involves py instead of ^a, and is ac- 
cordingly easier to work with in an Eulerian formulation 
The calculation of (62 1 is cumbersome if we use the 'di- 

because its non- 



rect' WKBJ solution for py, given by ( |52[ | 
trivially phase shifted with respect to the WKBJ solutions 
for p and px, given by (57a I and (57b I, making the tempo- 



ral average over one oscillation period somewhat ill deflned. 
However, if we exploit PPV conservation and express py in 
terms of p and px, see ( |58[ ), which is also found to be more 
accurate for e ~ 1 (see Fig.[4|, the calculation is trivial and 
we obtain 



(64) 



In Fig. [5] we show the angular momentum flux (F^) 
associated with the wave parts of solutions plotted in Fig. |4] 
In order to obtain the wave part of the numerical solution 
we have simply subtracted the balanced solution. 
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exp ( — 7r/4e) 





1/256 1/64 



1/16 1/4 1 

kyH/2TT 



Figure 5. Wave action for the wave sliown in Fig. |4] The black 
and red hne are the wave action for the wave part of the numerical 
and WKBJ solution, respectively. The dashed line indicates the 
large-T average value K3t. 



As expected, the WKBJ solutions agrees with the nu- 
merical solution remarkably well. Due to interference be- 
tween the forward and the backwards travelling waves actu- 
ally obtained in our rmmerical solution, the quantity defined 



by (621 oscillates. However, the average over one oscillation 



period approaches a constant non-zero value as is expected 
for a linear wave. 



5 DISCUSSION 

In this paper we have developed a theory of SD wave exci- 
tation in a rotating shear flow with turbulence which may 
result from the magneto-rotational instability but under the 
assumption that the magnetic field is too weak to affect the 
form of the waves significantly. We considered the commonly 
adopted shearing box model for which the fiow is subject to 
the boundary condition of periodicity in shearing coordi- 
nates (Goldr eich fc Lynden-Bell|1965 1. 

The main feature resulting from the shear is that wave 
excitation occurs through a sequence of regularly spaced 
swings as the wave changes from leading to trailing form. 



For a fixed azimuthal wave number ky, and a Keplerian ro- 
tation profile, the swings are separated by a time interval 
Sts — Toib/iqkyLx), where the orbital period Torb = 2tt/Q. 
For the optimal azimuthal wave number ky — k/c (see dis- 
cussion below) and Lx = H as an estimated radial correla- 
tion length of the turbulence, which should also be the min- 
imum box size required to capture its essential properties, 

it follows that Sts ~ Torb. 

The wave equations governing the excitation during a 
particular swing were found to depend on time alone and un- 
der the assumption that the important source terms causing 
the wave excitation are associated with the pseudo poten- 
tial vorticity, they could be solved to find the asymptotic 
wave form and net positive wave action produced. The form 
of the wave equations necessitated a WKBJ analysis in the 
complex plane. In this respect the formalism differs from 
shearing box analyses that adopt rigid or free boundary con- 
ditions, or which assume strictly harmonic forcing with ra- 
dial boundaries extended to infinity, rather than periodicity 
in shearing coordinates. In the former cases one can sepa- 



Figure 6. Exponential dependence of the WKBJ amplitudes on 
the azimuthal wave number ky for a Keplerian rotation law. 



rate out a harmonic time dependence and solve a problem 
in space for the wave amplitude (e.g. [Narayan et al. 11987] ). 

The analysis of the wave excitation process driven by 
pseudo potential vorticity carried out in this paper has sim- 
ilarities to an analysis of inertia-gravity waves excited in 

£2004) 



who 
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perform an analogous WKBJ analysis in the complex plane. 

The excitation process produces waves of equal ampli- 
tude propagating in opposite directions. As these waves are 
both trailing, by symmetry each produces an equal out- 
ward angular momentum flow. Even when the excitation 
process is linear, as the waves propagate away, the radial 
wave length shortens until shock dissipation eventually oc- 
curs (e.g. [Goodman fc Rafikov|2001| ). Thus waves are always 
likely to be seen to manifest nonlinear effects as the charac- 
teristic radial wave length shortens. When waves behave lin- 
early during the initial excitation, but subsequently undergo 
significant but not complete dissipation between successive 
swings, the rate of angular momentum transport can be es- 
timated as the wave action produced in single swing given 



by equation (63). This situation is found to occur in the 



simulations presented in paper II. 

An important parameter is the value of the az- 
imuthal/horizontal wave number, ky, for which the wave 
excitation is most favoured, which we define as being the 
value for which the wave amplitude produced is maximal. 



According to ( 52 1 and ( 55 1 the wave amplitude is a prod- 



uct of a known function of ky and the square of the Fourier 
amplitude of PPV at the time of swing. The latter quantity, 
being determined by the nonlinear hydromagnetic turbu- 
lence, cannot be found from the wave excitation calculations 
performed here. This aspect is discussed in paper II where 
relevant numerical simulations are performed and analysed. 
Here we shall anticipate results and assume that the PPV 
spectrum is relatively flat at small ky with the consequence 
that it does not affect the dependence of the wave amplitude 
on ky significantly. 

The wave amplitude depends on the azimuthal wave 
number ky through the parameter 

qilkyC 
kyC^ + K? 

in such a way that it is exponentially small for e <C 1, see 
(|53|) and (1561 . Given the fact that e is small both in the 
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small and the long azimuthal wave number limit, we deduce 
that wave excitation will be most effective near the optimal 
wave number 

for which e takes its maximum value 

tmax = qQ,/(2li). 

For a Keplerian disc with q = 3/2 and k = 1 we have 
fc^* = l/H and thus e^^^ = 3/4. We recaU that WKBJ 
theory gives very accurate results for values of e as large as 
this. For illustrative purposes we plot the exponential de- 
pendence of the wave amplitudes as a function of azimuthal 
wave number in Fig. [5] We see that amplitude of the excited 
wave falls off rapidly away from the optimal wave number. 

The arguments given above suggest that SD wave pro- 
duction will be most effective for ky ~ fcy^'- This is the 
longest possible azimuthal wave length for a box with 
Ly = 2nH as is commonly adopted. For boxes of this size 
and smaller wave production is expected to be most effec- 
tive at the longest azimuthal wave length. On the other hand 
once Ly exceeds 2nH the longest wave length is expected to 
no longer be the most effective. This is fully supported by 
the simulation results presented in paper II. In this paper 
we confirm the main features of the excitation process de- 
scribed here and verify the dominance of the pseudo poten- 
tial vorticity related source terms. Although the waves are 
observed to become nonlinear very soon after the initial ex- 
citation, the main features of the analysis presented here are 
confirmed. This suggests that useful extensions can be made 
to the analysis of wave excitation under more general condi- 
tions such as those that incorporate significant self-gravity. 
We plan to undertake these in the near future. 
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